% This script starts estimation of the Macro finance model using the SR approach
% By Martin M. Andreasen, August 2019
addpath(genpath(pwd))
close all
clear 
%clc
%% User settings
dataOption        = 1;          %1 for using ygap = UGAP,
                                %2 for using ygap = potential output in GDP
                                %3 for using ygap = Cooper and Priestley gap
startYear         = 1961;       %startYear for the estimation. Data available from 1961
modelNum          = 1;          %1 for QTSM without interest rate smoothing
                                %2 for QTSM with interest rate smoothing

% For inference
theta12OptimalOn  = 0;          %1 for setting LAMBDA in step 3 in an optimal manner. 0 for using the estimates from step 2
AVarTheta1On      = 0;          %1 for computing AVar for theta1 at step 1 and theta11 at step 3, else 0

% The optimizer:
optimizerStep1     = 3*0;         % optim = 0 for no optimization
optimizerStep3     = 3*0;         % optim = 1 for the CMAES
                                % optim = 2 for a gradient based optimizer (the LM version)
                                % optim = 3 for the Nelder-Mead.
                                % optim = 4, for first a gradient optimizer and the the Nelder-Mead
optim.numOptim    = 1;          % Number of optimizations - if optimizer = 3 or 4.
optim.sigma       = 0.10;       % step length for CMAES
optim.sigmaMax    = 0.20;       % max length for CMAES
optim.TolFun      = 1e-5;
optim.TolX        = 1e-5;
optim.MaxIter     = 4000;
optim.MaxFunEvals = 4000;

epsValue          = 1D-5;       % Stepsize for the standard errors of Theta1
AVarAdjThetaOn    = 0;          % 1 to adjust factor variance for Theta1, else 0

h0RegimeOn        = 0;          %1 for switching in h0 under P, else 0
hxRegimeOn        = 1;          %1 for switching in hx under P, else 0

%% Default setting
numCPUs           = feature('numcores'); % Number of CPUs
fullPHIon         = 1;          %1 for a full phi, else 0 for diagonal phi
nm                = 2;          %Number of macro factors
nn                = 0;          %Number of latent factors factors
loadDefaultSetting

% Loading the data:
loadData = getData(dataOption,startYear);
%% Starting values for the estimated parameters
theta1 = initializeTheta1(fullPHIon,nm,nn);

% Impose restrictions on phi from ReStud paper
theta1 = rmfield(theta1,'phi14');
theta1 = rmfield(theta1,'phi23');
theta1 = rmfield(theta1,'phi32');
theta1 = rmfield(theta1,'phi34');
theta1 = rmfield(theta1,'phi41');
theta1 = rmfield(theta1,'phi43');

if modelNum == 1  %
    theta1 = rmfield(theta1,'rhor');
    theta1ValuesStep1 = [   0.000007697957    0.088798954825    0.000001448479    0.006434745357    0.025271178113  ...
        0.003801472813    0.004258045841   -0.485207099637   -0.491153447384   -0.021787206525  ...
        -0.000928735436    0.002855108937   -0.024156823397    0.000103525349    0.000830331327  ...
        0.023919846878   -0.087674048434    0.049990615980    0.018135323127   -0.000302856709 ...
        0.039970173130    0.002512837718   -0.004243300828    0.002117672196    0.000003903177  ]';
    % Q = 0.069862979582060, done by CMAES
    
    % For step 3
    if h0RegimeOn == 0 && hxRegimeOn == 0
        theta1ValuesStep3 = [   -0.000054267465    0.062141749544    0.000001029203    0.020190705171    0.023033327077  ...
            0.002630766803    0.004823063396   -0.252093024059   -0.361085717407    0.015628416992  ...
            -0.000725748632    0.002047542392   -0.015202204828    0.004076399520    0.000672800837  ...
            0.001655995736    0.000039057451    0.001943359316   -0.006279208110   -0.001922424891  ...
            0.011801466352   -0.000719350910   -0.000207358753    0.000926326930    0.000712422918 ]';
        % Q =   0.073725835300225, done by CMAES
    elseif h0RegimeOn == 1 && hxRegimeOn == 0
        theta1ValuesStep3 = [ -0.000052747743    0.062111513491    0.000001081952    0.020199920308    0.023027515066  ...
            0.002629878536    0.004822195881   -0.251959149979   -0.361127843815    0.015658391469  ...
            -0.000725639430    0.002047666229   -0.015212567435    0.004076203140    0.000672794230  ...
            0.001634277605    0.000031876403    0.001942920152   -0.006491309075   -0.001961645153 ...
            0.011653109936   -0.000728301798   -0.000209827337    0.000919178644    0.000711854713 ]';
        % Q = 0.073727124828977, done
    elseif h0RegimeOn == 0 && hxRegimeOn == 1
        theta1ValuesStep3 = [   -0.000052529151    0.061893644107    0.000001015031    0.020127398256    0.023015775710  ...
            0.002649581929    0.004818868026   -0.249851076249   -0.362188831135    0.015830591983  ...
            -0.000727219081    0.002044829664   -0.015189659630    0.004084626447    0.000671770332  ...
            0.001634667604    0.000033630677    0.001939529654   -0.006542854015   -0.001840507312 ...
            0.011472572799   -0.000733829336   -0.000202007946    0.000909291558    0.000711139421 ]';
        % Q= 0.073727371238604, done 
    elseif h0RegimeOn == 1 && hxRegimeOn == 1
        % Change in intercepts are not significant
        theta1ValuesStep3 = [     -0.000105253887    0.414352586444    0.494286583897    0.301285280938    0.000001093386  ...
            0.073390765022    0.287349341969   -0.041571247513   -0.003714475854    0.389062513331  ...
            -0.088959177209    0.037238279387    0.085480399328    0.026462625781   -0.015214482809  ...
            0.001693778664    0.000059919995    0.001949913807   -0.006211423343   -0.001786126762 ....
            0.011687749617   -0.000721208403   -0.000199725004    0.000921912262    0.000711221986  ]';
        % Q = 0.096928521607404, done by CMAES
    end
elseif modelNum == 2
    % Step1
    theta1ValuesStep1 = [     0.000002147961    0.000053633462    0.088207267823    0.000002127190    0.006429962185  ...
        0.025056306897    0.003774613632    0.004390757339   -0.481685799885   -0.498985394823  ...
        -0.021947420144   -0.000889881215    0.002824329345   -0.023938927490   -0.000020696901  ...
        0.000810077481    0.023731220181   -0.086298933061    0.049861640544    0.017870718919  ...
        -0.000540883622    0.039536445323    0.002453193621   -0.004176233326    0.002074913661  ...
        0.000008088638]';
    % QStep1 = 0.069874703757766, done
    
    % For step 3
    if h0RegimeOn == 0 && hxRegimeOn == 0
        
        % Q =   
    elseif h0RegimeOn == 1 && hxRegimeOn == 0
        
    elseif h0RegimeOn == 0 && hxRegimeOn == 1
        theta1ValuesStep3 = [    -0.000006740469    0.000012251801    0.064388794372    0.000001029382    0.021418625542  ...
            0.023759336398    0.002685638941    0.004899239166   -0.256190028514   -0.338525350799  ...
            0.014477814377   -0.000758245029    0.002224218891   -0.015623441106    0.004127740071  ...
            0.000723405862    0.001634658460    0.000033401458    0.001939559263   -0.006544495627  ...
            -0.001852831827    0.011489231888   -0.000727163053   -0.000200410401    0.000909873435  ...
            0.000695994319 ]';
        % Q = 0.073840225866708,done
    elseif h0RegimeOn == 1 && hxRegimeOn == 1
        
    end
end
theta1           = values2struct(theta1ValuesStep1,fieldnames(theta1));
theta1StartStep3 = values2struct(theta1ValuesStep3,fieldnames(theta1));

% Getting bounds and Insigma
[upperBounds,lowerBounds,Insigma] = theta1BoundsInsigma(nm,nn);

% Defining the struct with the setup specification for step 1
constructSetupStep1

%% Computing the objective function theta1
tic
[Q,output,errorMes] = objectFunc(struc2values(theta1),setupStep1);
toc
%% Step 1 of the SR approach
% The optimization
optim.optimizer = optimizerStep1;
[QStep1,theta1Step1,outputStep1] = step1SR(theta1,setupStep1,optim);
theta1 = theta1Step1;

% The standard errors
AVarTheta1Step1  = getAVarTheta1andFactors(theta1Step1,setupStep1,epsValue);
%plotFactors(setupStep1,outputStep1,AVarTheta1Step1)
%% Step 2 of the SR approach
[theta2Step2,AVarTheta2Step2] = step2SR_threshold_macro(AVarTheta1Step1,outputStep1,bandwidthTstep2,setupStep1.econAct,h0RegimeOn,hxRegimeOn,numBootStep2);
%mprStep2 = getMPR(outputStep1,theta2Step2,AVarTheta2Step2,nm*2+nn,h0RegimeOn,hxRegimeOn);
%% Step 3 of the SR approach
namesTheta12 = cell((2*nm+nn)*(2*nm+nn+1)/2,1);
idx = 0;
for i=1:2*nm+nn
    for j=1:i
        idx = idx + 1;
        namesTheta12{idx} = ['sigma',num2str(i),num2str(j)];
    end
end
[theta12Step3,AVarTheta12Step3] = step3SR_theta12(namesTheta12,theta1Step1,AVarTheta1Step1,theta2Step2,AVarTheta2Step2,theta12OptimalOn);

% Estimation of theta11 - using Robust version
if exist('theta1StartStep3','var')
    theta1Start = theta1StartStep3;
    if AVarTheta1On  == 1
        disp('for standard errors')
        theta1Start = rmfield(theta1Start,'phi22');
        setupStep1.calibrateTheta1.phi22 = theta1StartStep3.phi22;
        setupStep1.selectTheta1          = fieldnames(theta1Start);       
        setupStep1.upperBoundsValues    = struc2values(upperBounds,setupStep1.selectTheta1);
        setupStep1.lowerBoundsValues    = struc2values(lowerBounds,setupStep1.selectTheta1);
        setupStep1.InsigmaValues        = struc2values(Insigma,setupStep1.selectTheta1);
     end
else
    theta1Start = theta1Step1;
end
optim.optimizer = optimizerStep3;
[QStep3,theta11Step3,outputStep3,setupStep3] = step3SR_theta11(theta12Step3.Robust,namesTheta12,theta1Start,...
    setupStep1,upperBounds,lowerBounds,Insigma,optim);

% Standard errors of theta11: The asymptotic standard errors 
AVarTheta11Step3 = getAVarTheta11Step3(namesTheta12,theta11Step3,AVarTheta1Step1,AVarTheta2Step2,AVarTheta12Step3,...
    setupStep1,setupStep3,epsValue/10,theta12OptimalOn);

% Re-estimating P dynamics
[theta2Step3,AVarTheta2Step3] =step2SR_threshold_macro(AVarTheta11Step3,outputStep3,bandwidthTstep2,setupStep3.econAct,h0RegimeOn,hxRegimeOn,numBootStep2);

% Enforce consistency of sigma between Q and P
posTheta12Step2 = zeros(1,length(namesTheta12));
for i=1:length(namesTheta12)
    theta2Step3.(namesTheta12{i}) = theta2Step2.(namesTheta12{i});
    posTheta12Step2(1,i) = find(strcmp(AVarTheta2Step3.names,namesTheta12(i)));
end
AVarTheta2Step3.VarRobust(posTheta12Step2,posTheta12Step2) = AVarTheta2Step2.VarRobust(posTheta12Step2,posTheta12Step2); 

mprStep3 = getMPR(outputStep3,theta2Step3,AVarTheta2Step3,nm*2+nn,h0RegimeOn,hxRegimeOn);
testRegimeStep3 = waldTestRegimeShift(theta2Step3,AVarTheta2Step3,2*nm+nn);

%% Do the yield curve decomposition
if h0RegimeOn == 0 && hxRegimeOn == 0
    yieldTP = yieldCurveDecom(setupStep3,outputStep3,theta2Step3);
else
    numSimTP = 1D4;
    yieldTP  = yieldCurveDecom_threshold(setupStep3,outputStep3,theta2Step3,numSimTP);
end

% mean(yieldTP.termPremia(end,:))
% std(yieldTP.termPremia(end,:))
%% Simulating the model
thresholdLevel          = AVarTheta2Step3.thresholdLevel;
momentsData             = getMomentsData(loadData,outputStep3,thresholdLevel);
momentsModel            = simulator_threshold(theta2Step3,outputStep3,thresholdLevel,numSim);
momentsModel_finite     = simulator_threshold_finite(theta2Step3,outputStep3,thresholdLevel,1000);
expectedExcessReturn.h1 = getExpectedExcessReturns(loadData,setupStep3,outputStep3,theta2Step3,1);
expectedExcessReturn.h3 = getExpectedExcessReturns(loadData,setupStep3,outputStep3,theta2Step3,3);

%% Do plots
plotFactors(setupStep3,outputStep3,AVarTheta11Step3)
plotFit(setupStep3,outputStep3)
plotUnconditionalMoments(momentsModel,momentsData)
plotCSregressions(momentsModel,momentsData)
plotForecastRegressions(momentsModel,momentsData)
plotYieldCurveDecom(yieldTP,setupStep3)

%% Displaying the results
outTheta1          = displayResults(theta1Step1,AVarTheta1Step1.VarRobust);
outTheta11         = displayResults(theta11Step3,AVarTheta11Step3.VarRobust);
theta2Step3Reduced = reduceTheta2(theta2Step3,AVarTheta2Step3.names);
outTheta2Reduced   = displayResults(theta2Step3Reduced,AVarTheta2Step3.VarRobust);

disp('optimum at step 1')
printToScreen(struct2array(theta1Step1));
disp('optimum at step 3')
theta1Step3Start  = theta11Step3;
for i = 1:size(namesTheta12,1)
    theta1Step3Start.(namesTheta12{i}) = setupStep3.calibrateTheta1.(namesTheta12{i});
end
printToScreen(struct2array(theta1Step3Start));

%% Saving the results in a mat file
if modelNum == 1
    modelName = 'QTSM_nn0_ReStudPhiQ_';
elseif modelNum == 2
    modelName = 'QTSMr_nn0_ReStudPhiQ_';
end
filename = [modelName,num2str(startYear),'_SE',num2str(epsValue),'_h0RegimeOn',num2str(h0RegimeOn),'_','hxRegimeOn',num2str(hxRegimeOn),...
    '_numBootStep2_',num2str(numBootStep2)];
save([pwd,'\Results\',filename,'.mat'],'QStep1','theta1Step1','outputStep1','AVarTheta1Step1',...
     'theta2Step2','AVarTheta2Step2',...
     'theta11Step3','outputStep3','setupStep3','AVarTheta11Step3','theta2Step3','AVarTheta2Step3',...
     'yieldTP','momentsModel','momentsModel_finite','momentsData','mprStep3','testRegimeStep3','expectedExcessReturn')



